function Fisher = ecmmvnrfish(Data, Design, Covar, Method, MatrixFormat, CovarFormat)
%ECMMVNRFISH Fisher information for multivariate normal regression model.
%	Fisher information matrix based on current maximum likelihood or least-squares parameter
%	estimates that account for missing data.
%
%		Fisher = ecmmvnrfish(Data, Design, Covar);
%		Fisher = ecmmvnrfish(Data, Design, Covar, Method, MatrixFormat, CovarFormat);
%
% Inputs:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a NUMSERIES-dimensional random
%		vector. Missing values are represented as NaNs. Only samples that are entirely NaNs are
%		ignored (to ignore samples with at least one NaN, use mvnrfish).
%	Design - Either a matrix or a cell-array to handle two distinct model structures. First, if
%		NUMSERIES = 1, Design can be a NUMSAMPLES x NUMPARAMS matrix with known values. This is the
%		"standard" form for regression on a single data series. Alternatively, for any
%		NUMSERIES >= 1, Design can be a cell array of either 1 or NUMSAMPLES cells, where each cell
%		contains a NUMSERIES x NUMPARAMS matrix of known values. If Design has a single cell, then
%		it is assumed to be the same Design matrix for each sample. Otherwise, Design must contain
%		individual Design matrices for each sample.
%	Covar - NUMSERIES x NUMSERIES matrix of estimates for the covariance of the residuals of the
%		regression.
%
% Optional Inputs:
%	Method - String to identify method of calculation for information matrix. The choices are:
%		'hessian' - Default method. Use the expected hessian of the observed log-likelihood
%			function. This method is recommended since the resultant standard errors incorporate the
%			increased uncertainties due to missing data.
%		'fisher' - Use the fisher information matrix.
%	MatrixFormat - String that identifies parameters to be included in the Fisher information
%		matrix. The choices are:
%		'full' - Default format. Compute the full Fisher information matrix for both model and
%			covariance parameter estimates.
%		'paramonly' - Compute only components of the Fisher information matrix associated with the
%			model parameter estimates.
%	CovarFormat - String that specifies the format for the input covariance matrix. The choices are:
%		'full' - Default method. Compute the full covariance matrix.
%		'diagonal' - Treat the covariance matrix as a diagonal matrix.
%
% Outputs:
%   Fisher - TOTALPARAMS x TOTALPARAMS Fisher information matrix or Hessian matrix based on current
%		parameter estimates, where
%			TOTALPARAMS = NUMPARAMS + NUMSERIES * (NUMSERIES + 1)/2
%		if MatrixFormat = 'full' and
%			TOTALPARAMS = NUMPARAMS
%		if MatrixFormat = 'paramonly'.
%
% WARNING: If calculating the full Fisher information or Hessian matrix, this function is VERY slow.
%
% See also ECMMVNRSTD, ECMMVNRMLE.

%	Copyright 2005-2007 The MathWorks, Inc.
%	$Revision: 1.1.6.3 $ $Date: 2007/05/10 13:44:56 $

% Step 1 - check arguments

if nargin < 6 || isempty(CovarFormat)
	CovarFormat = 'FULL';
end
if nargin < 5 || isempty(MatrixFormat)
	MatrixFormat = 'FULL';
end
if nargin < 4 || isempty(Method)
	Method = 'HESSIAN';
end
if nargin < 3
	error('Finance:ecmmvnrfish:MissingInputArg', ...
		'Missing required input arguments Data, Design, or Covar.');
end
if isempty(Data)
	error('Finance:ecmmvnrfish:EmptyDataArray', ...
		'The required input argument Data is empty.');
end
if isempty(Design)
	error('Finance:ecmmvnrfish:EmptyDesignArray', ...
		'The required input argument Design is empty.');
end
if isempty(Covar)
	error('Finance:ecmmvnrfish:EmptyCovar', ...
		'The required input argument Covar is empty.');
end

%[NumSamples, NumSeries, NumParams] = ...
%	checkmvnrsetup(Data, Design, [], Covar);

[NumSamples, NumSeries] = size(Data);
if iscell(Design)
	if (numel(Design) == 1)
		SingleDesign = true;
	else
		SingleDesign = false;
	end
	NumParams = size(Design{1},2);
else
	SingleDesign = false;
	NumParams = size(Design,2);
end

if ~all(size(Covar) == [NumSeries, NumSeries])
	error('Finance:ecmmvnrfish:InconsistentDims', ...
		'The covariance matrix Covar has wrong dimensions.');
else
	[CholCovar, CholState] = chol(Covar);
	if CholState > 0
		error('Finance:ecmmvnrfish:NonPosDefCov', ...
			'Covariance matrix is not positive-definite.');
	end
end

% Step 2 - initialization

if ~any(strcmpi(Method,{'FISHER','HESSIAN'}))
	warning('Finance:ecmmvnrfish:UnknownMethodString', ...
		'Unknown Method string. Will use default HESSIAN.');
	Method = 'HESSIAN';
end

if ~any(strcmpi(MatrixFormat,{'PARAMONLY','FULL'}))
	warning('Finance:ecmmvnrfish:UnknownFormatString', ...
		'Unknown MatrixFormat string. Will use default FULL.');
	MatrixFormat = 'FULL';
end

if strcmpi(MatrixFormat,'PARAMONLY')
	TotalParams = NumParams;
else
	if strcmpi(CovarFormat,'FULL')
		TotalParams = NumParams + (NumSeries * (NumSeries + 1))/2;
	else
		TotalParams = NumParams + NumSeries;
	end
end

Fisher = zeros(TotalParams,TotalParams);

% Step 3 - calculate fisher information matrix (not hessian)

if strcmpi(Method,'FISHER')
	
	% Step 4 - do partials wrt Mean

	Count = 0;
	TestMatrix = zeros(NumParams,NumParams);
	if iscell(Design)
		if SingleDesign
			A = CholCovar' \ Design{1};
			for k = 1:NumSamples
				TestMatrix = TestMatrix + A'*A;
			end
			Count = NumSamples;
		else
			for k = 1:NumSamples
				if ~all(isnan(Data(k,:)))
					Count = Count + 1;
					A = CholCovar' \ Design{k};
					TestMatrix = TestMatrix + A'*A;
				end
			end
		end
	else
		for k = 1:NumSamples
			if ~all(isnan(Data(k,:)))
				Count = Count + 1;
				A = CholCovar' \ Design(k,:);
				TestMatrix = TestMatrix + A'*A;
			end
		end
	end
	TestMatrix = (1.0/Count) .* TestMatrix;
	Fisher(1:NumParams,1:NumParams) = TestMatrix;

	% Step 5 - do partials wrt Covar

	if strcmpi(MatrixFormat,'FULL')
		if strcmpi(CovarFormat,'FULL')
			InvCovar = inv(Covar);

			GradC1 = zeros(NumSeries,NumSeries);
			GradC2 = zeros(NumSeries,NumSeries);

			i = NumParams;
			for i1 = 1:NumSeries
				for j1 = 1:i1
					i = i + 1;

					GradC1(i1,j1) = 1;						% do dC/dtheta(i)
					GradC1(j1,i1) = 1;

					j = NumParams;
					for i2 = 1:NumSeries
						for j2 = 1:i2
							j = j + 1;

							if (j <= i)
								GradC2(i2,j2) = 1;			% do dC/dtheta(j)
								GradC2(j2,i2) = 1;

								Temp1 = InvCovar*GradC1;
								Temp2 = InvCovar*GradC2;

								Fisher(i,j) = 0.5*trace(Temp1*Temp2);
								Fisher(j,i) = Fisher(i,j);                        

								GradC2(i2,j2) = 0;			% undo dC/dtheta(j)
								GradC2(j2,i2) = 0;
							end
						end
					end

					GradC1(i1,j1) = 0.0;                    % undo dC/dtheta(i)
					GradC1(j1,i1) = 0.0;
				end
			end
		else
			CDiag = 1 ./ diag(Covar);
			FDiag = 0.5 * CDiag .^2;
			Fisher(NumParams*TotalParams + NumParams+1:TotalParams+1:end) = FDiag;
		end
	end
	
% Step 6 - calculate hessian (not fisher information matrix)

else
	% Step 7 - main loop over data records

	if strcmpi(CovarFormat,'FULL')
		InvCovar = inv(Covar);
	else
		InvCovar = diag(1 ./ diag(Covar));
	end
	
	Map = zeros(NumSeries,1);
	Count = 0;

	for kk = 1:NumSamples

		% Step 8 - determine and map available data in current record

		Map(:) = 0;
		Available = 0;
		for ii = 1:NumSeries
			if isnan(Data(kk,ii))
				Map(ii) = 0;
			else
				Map(ii) = 1;
				Available = Available + 1;
			end
		end

		if Available > 0                        % skip over empty records
			Count = Count + 1;

			% Step 9 - construct covariance matrix subarrays

			if iscell(Design)
				if SingleDesign
					SubDesign = Design{1};
				else
					SubDesign = Design{kk};
				end
			else
				SubDesign = Design(kk,:);
			end
			SubCovar = Covar;

			if Available < NumSeries
				for ii = NumSeries:-1:1
					if Map(ii) == 0
						SubDesign(ii,:) = [];
						SubCovar(:,ii) = [];
						SubCovar(ii,:) = [];
					end
				end
				if strcmpi(CovarFormat,'full')
					InvSubCovar = inv(SubCovar);
				else
					InvSubCovar = diag(1 ./ diag(SubCovar));
				end
			else
				InvSubCovar = InvCovar;
			end

			% Step 10 - do partials wrt Mean for current data record

			TempMatrix = SubDesign' * InvSubCovar * SubDesign;
			Fisher(1:NumParams,1:NumParams) = Fisher(1:NumParams,1:NumParams) + TempMatrix;
			
			% Step 11 - do partials wrt Covar for current data record

			if strcmpi(MatrixFormat,'FULL')
				if strcmpi(CovarFormat,'FULL')
					GradC1 = zeros(Available,Available);
					GradC2 = zeros(Available,Available);

					p1 = 0;
					i = NumParams;
					for i1 = 1:NumSeries
						if Map(i1) > 0
							p1 = p1 + 1;
						end

						q1 = 0;
						for j1 = 1:i1
							i = i + 1;

							if Map(j1) > 0
								q1 = q1 + 1;
							end

							if (Map(i1) > 0) && (Map(j1) > 0)
								GradC1(p1,q1) = 1;
								GradC1(q1,p1) = 1;
							end

							p2 = 0;
							j = NumParams;
							for i2 = 1:NumSeries
								if Map(i2) > 0
									p2 = p2 + 1;
								end

								q2 = 0;
								for j2 = 1:i2
									j = j + 1;

									if Map(j2) > 0
										q2 = q2 + 1;
									end

									% dC/dtheta(i) = dC/dC(i1,j1)
									% dC/dtheta(j) = dC/dC(i2,j2)

									if (j <= i) && (Map(i1) > 0) && (Map(j1) > 0) && (Map(i2) > 0) && (Map(j2) > 0)
										GradC2(p2,q2) = 1;
										GradC2(q2,p2) = 1;

										Temp1 = InvSubCovar*GradC1;
										Temp2 = InvSubCovar*GradC2;

										Fisher(i,j) = Fisher(i,j) + 0.5*trace(Temp1*Temp2);
										Fisher(j,i) = Fisher(i,j);

										GradC2(p2,q2) = 0;
										GradC2(q2,p2) = 0;
									end
								end
							end

							if (Map(i1) > 0) && (Map(j1) > 0)   % undo dC/dtheta(i)
								GradC1(p1,q1) = 0;
								GradC1(q1,p1) = 0;
							end
						end
					end
				else
					GradC1 = zeros(Available,Available);

					p1 = 0;
					i = NumParams;
					for i1 = 1:NumSeries
						i = i + 1;

						if Map(i1) > 0
							p1 = p1 + 1;

							GradC1(p1,p1) = 1;
						
							Temp1 = InvSubCovar*GradC1;

							Fisher(i,i) = Fisher(i,i) + 0.5*trace(Temp1*Temp1);
							
							GradC1(p1,p1) = 0;
						end
					end
				end
			end
		end
	end

	% Step 12 - normalize hessian

	Fisher = (1.0/Count) .* Fisher;
end
